function mu=equilibrium(nw,nm,Phira)

C=nm*nw+nm+nw; %number of unknowns

%% Equality constraints Aeq*mu=beq
%(1) Men: sum of mu across women (including singles)=1
coordrowM=kron(1:1:nm, ones(1,nw+1)); %1x(nm*(nw+1))
temp=zeros(nw+1, nm);
vec = (1:nm:(nm*nw+1)).';
for ii=1:nm
    temp(:,ii) = vec;
    vec = vec + 1;
end
coordcolumnM=temp(:).'; %1x(nm*(nw+1))
fillAeqM=ones(1,nm*(nw+1)); %1x(nm*(nw+1))

%(2) Women: sum of mu across men (including singles)=1
coordrowW=kron(coordrowM(end)+1:1:coordrowM(end)+nw, ones(1,nm+1)); %1x(nw*(nm+1))
temp=zeros(nm+1,nw);
vec = (1:nm).';
for jj = 1:nw
    temp(1:end-1,jj) = vec;
    vec = vec + nm;
end
temp(end,:) = nm*nw+nm + (1:nw); 
coordcolumnW=temp(:).'; %1x(nw*(nm+1))
fillAeqW=ones(1,nw*(nm+1)); %1x(nw*(nm+1))

%(3) Compose
coordrow=[coordrowM coordrowW]; 
coordcolumn=[coordcolumnM coordcolumnW];
fillAeq=[fillAeqM fillAeqW];
Aeq=sparse(coordrow, coordcolumn, fillAeq, nm+nw, C); %(nm+nw)xC
beq=ones(nm+nw,1); %(nm+nw)x1


%% Objective function
f=-Phira;  %Cx1

%% MOSEK (Slower and does not exactly satisfy the integer constraints)
%prob.blx=zeros(C,1); %lower bound unknowns
%prob.ulx=ones(C,1); %upper bound unknowns
%prob.c=-f; %objective function
%param_MOSEK.MSK_IPAR_LOG = 0;
%prob.a=Aeq;
%prob.blc=beq;
%prob.buc=beq;
%[~,result]     = mosekopt('minimize echo(0)',prob, param_MOSEK);
%mu=result.sol.itr.xx;

%% Gurobi
model.A=Aeq;
model.obj=f;
model.sense=repmat('=', size(Aeq,1),1);
model.rhs=beq; 
model.ub=ones(C,1); 
model.lb=zeros(C,1); 
params.outputflag = 0; 
result_GUROBI=gurobi(model, params);
mu=result_GUROBI.x; %Cx1

end
